Simulating RCTs as Priors (Models)
Arthur M. Albuquerque
September 15, 2021
# Ensures the package "pacman" is installed
if (!require("pacman")) install.packages("pacman")
pacman::p_load(tidyverse, # data wrangling + plotting
gt, # to create tables
brms, # to fit Bayesian models
here, # reproducible file paths
rio, # to import files
ggdist, # to plot distributions
tidybayes, # to plot distributions
PNWColors,
flextable, # table
patchwork) # arrange plots
# Load data
d_logOR = readRDS(
here("03_updated_analyses_who", "output", "data", "effect_sizes.Rds")
)
# Load main model
m1 = readRDS(here("03_updated_analyses_who", "output", "fits", "m1.Rds"))
# Load function
source(here("03_updated_analyses_who/functions/diag_plot.R"))In this file, we will discuss the rationale and fit models in which we use simulated randomized controlled trials as priors.
Details
The main Bayesian meta-regression random-effect model is:
\[ \begin{align*} y_i & \sim Normal(\theta_i, \sigma_i^2) \tag{Likelihood} \\ \theta_i & \sim Normal(\mu, \tau^2)\\ \mu & = \alpha + \beta_{SOO} x_i + \beta_{NIV}x_i\\ \\ \alpha & \sim \operatorname{Normal}(0, 0.82^2) \tag{Priors} \\ \beta_{SOO}, \beta_{NIV} & \sim \operatorname{Normal}(0, 1.5^2) \\ \tau & \sim \operatorname{Half-Normal}(0.5^2) \\ \end{align*} \]
where \(y_i\) is the observed mean log odds ratio of tocilizumab versus control and \(\sigma_i^2\) is the known sampling variance in study \(i\). \(\alpha\) is the intercept, which represents the overall effect of tocilizumab in patients on invasive mechanical ventilation. \(\beta_{SOO}\) is the difference between patients on simple oxygen only (SOO) and patients on invasive mechanical ventilation (\(\alpha\)). \(\beta_{NIV}\) is the difference between patients on noninvasive ventilation (NIV) and patients on invasive mechanical ventilation (\(\alpha\)). Both coefficients are multiplied by \(x\), which is dummy-coded. Lastly, \(\tau^2\) represents the between-study heterogeneity.
### Overall mean odds ratio in WHO's meta-analysis on tocilizumab in patients
# using corticosteroids:
# The mean odds ratio in patients using corticosteroids (toci vs control)
# was 0.77 (page 14 in Supplement 2 doi:10.1001/jama.2021.11330)
mean_OR_WHO = 0.77 Here, we will fit 6 models using different priors for \(\alpha\). Each prior will be tailored to simulate a randomized controlled trial (RCT) of tocilizumab versus standard of care (control). We want to simulate RCTs with the following total sample sizes: \(200, 500, 1000, 1500, 2000, 4000\).
All priors will be normally distributed and will have a mean of -0.26. This value is the equal to log(0.77), which was chosen based on WHO’s meta-analysis (page 14 in Supplement 2). This is the mean odds ratio of tocilizumab vs. control in patients using corticosteroids (overall results).
To calculate the standard deviation of each prior, we need to define certain aspects of these “simulated randomized controlled trials”:
control_risk = 0.43
# Now, we need to calculate the risk in the tocilizumab arm based on the odds ratio
# and control risk mentioned above
# Reference for formula: Box 1 in https://doi.org/10.1016/j.jclinepi.2020.08.019
toci_risk = (control_risk*mean_OR_WHO)/(control_risk*(mean_OR_WHO - 1) + 1)We assume equal allocation in both treatment arms.
Adapting from the suggestions in the GRADE guidelines, we found a striking discrepancy between the control mortality risk in this data (52%) in comparison to another previously published meta-analysis (34% in patients on IVM and using corticosteroids). Thus, we have decided to use 43% (arithmetic mean between 34 and 52) as our reference risk in the IVM subgroup..
The mortality risk in the tocilizumab was calculated using the following formula (Box 1 in Doi et al., 2020):
\[R_{T} = \frac{R_{C} O R}{R_{C}\left(O R-1\right)+1}\]
where \(R_{T}\) is the mortality risk in the tocilizumab group, \(R_{C}\) is the mortality risk in the control group and \(OR\) is the odds ratio mentioned above. Thus, the tocilizumab risk is equal to 0.37.
Thus, we are simulating RCTs with mean OR equal to 0.77, control risk mortality of 0.43, and tocilizumab risk of 0.37.
Based on these values, we can estimate the standard deviation (\(SD\)) with the following formula:
\[SD=\sqrt{\frac{1}{a+\frac{1}{2}}+\frac{1}{b+\frac{1}{2}}+\frac{1}{c+\frac{1}{2}}+\frac{1}{d+\frac{1}{2}}}\]
where \(a\), \(b\), \(c\) and \(d\) are number of events and follow this 2x2 table:
tribble(
~"Event", ~"Tocilizumab", ~"Control",
"Death", "a", "c",
"No death", "b", "d"
) %>%
flextable() %>%
autofit()Event | Tocilizumab | Control |
Death | a | c |
No death | b | d |
As similarly shown in the supplementary material of Higgins and Spiegelhalter, 2002, we can estimate these values as:
- \(a = R_{T}SS_{T}\)
- \(b = SS_{T} - R_{T}SS_{T}\)
- \(c = R_{C}SS_{C}\)
- \(d = SS_{C} - R_{C}SS_{C}\)
where \(SS_{T}\) and \(SS_{C}\) are the sample sizes in the tocilizumab and control arms, respectively. As mentioned above, we assume equal allocation in both treatment arms, thus \(SS_{T} = SS_{C}\).
# Function to calculate the standard deviation of odds ratio from simulated RCT ,
# based on risk and total sample size in each treatment arm
# Here we assume both arms have the same number of total patients (total_arm)
# Adapted from:
# Appendix in https://doi.org/10.1093/ije/31.1.96
# +
# Appendix in https://doi.org/10.1016/j.jclinepi.2008.07.006
SD_fun = function(risk_toci, # Risk in tocilizumab arm
risk_control, # Risk in control arm
total_arm){ # Total sample size in each arm
# Tocilizumab
a = risk_toci*total_arm # Deaths
b = total_arm - risk_toci*total_arm # No deaths
# Control
c = risk_control*total_arm # Deaths
d = total_arm - risk_control*total_arm # No deaths
# Standard deviation =
sqrt(1/(a + 1/2) + 1/(b + 1/2) + 1/(c + 1/2) + 1/(d + 1/2))
}# Lastly, we should apply the function above to calculate the standard deviation
# of multiple simulated RCTs will the same risk_toci and risk_control, but
# different amounts of patients
df =
tibble(risk_toci = toci_risk,
risk_control = control_risk,
total_arm = c(100, 250, 500, 750, 1000, 2000)) %>%
mutate(SD = SD_fun(risk_toci, risk_control, total_arm))Finally, we can estimate the \(SD\) based on the 6 different sample sizes mentioned above:
df %>%
summarise("Sample size in each treatment arm" = total_arm,
"Total sample size" = 2*total_arm,
" SD" = round(SD,2)) %>%
flextable() %>%
autofit()Sample size in each treatment arm | Total sample size | SD |
100 | 200 | 0.29 |
250 | 500 | 0.18 |
500 | 1,000 | 0.13 |
750 | 1,500 | 0.11 |
1,000 | 2,000 | 0.09 |
2,000 | 4,000 | 0.06 |
In summary, there are 4 parameters in these models, which are \(\alpha\), \(\beta_{SOO}\), \(\beta_{NIV}\), and \(\tau\). The priors for the latter three parameters will be the same in every model:
\[ \begin{align*} \beta_{SOO}, \beta_{NIV} & \sim \operatorname{Normal}(0, 1.5^2)\\ \tau & \sim \operatorname{Half-Normal}(0.5^2) \end{align*} \]
On the other hand, the priors for \(\alpha\) are described as:
\[ \begin{align*} \alpha & \sim \operatorname{Normal}(\mu, SD^2) \end{align*} \]
where \(\mu\) is equal to -0.26 in every model, but \(SD\) range from 0.29 to 0.06, as described in the table above.
Let’s fit these 6 new models:
# Priors
intercept100 = prior_string(
paste0(
"normal(",
log(mean_OR_WHO), # Mean, log(0.77)
",",
df[[1,4]], # SD, total_arm = 100
")"),
class = "b", coef = "Intercept"
)
# b_Intercept ~ normal(-0.261364764134408,0.28303369902181)
intercept250 = prior_string(
paste0(
"normal(",
log(mean_OR_WHO), # Mean, log(0.77)
",",
df[[2,4]], # SD, total_arm = 250
")"),
class = "b", coef = "Intercept")
# b_Intercept ~ normal(-0.261364764134408,0.179552672191073)
intercept500 = prior_string(paste0(
"normal(",
log(mean_OR_WHO), # Mean, log(0.77)
",",
df[[3,4]], # SD, total_arm = 500
")"),
class = "b", coef = "Intercept")
# b_Intercept ~ normal(-0.261364764134408,0.127092510741972)
intercept750 = prior_string(paste0(
"normal(",
log(mean_OR_WHO), # Mean, log(0.77)
",",
df[[4,4]], # SD, total_arm = 750
")"),
class = "b", coef = "Intercept")
# b_Intercept ~ normal(-0.261364764134408,0.103805945961058)
intercept1000 = prior_string(paste0(
"normal(",
log(mean_OR_WHO), # Mean, log(0.77)
",",
df[[5,4]], # SD, total_arm = 1000
")"),
class = "b", coef = "Intercept")
#b_Intercept ~ normal(-0.261364764134408,0.0899139032236548)
intercept2000 = prior_string(paste0(
"normal(",
log(mean_OR_WHO), # Mean, log(0.77)
",",
df[[6,4]], # SD, total_arm = 2000
")"),
class = "b", coef = "Intercept")
# b_Intercept ~ normal(-0.261364764134408,0.0635949873256457)
SOO = prior(normal(0, 1.5), class = "b", coef = "oxygenlow") # beta_Simple oxygen only
NIV = prior(normal(0, 1.5), class = "b", coef = "oxygenNIV") # beta_NIV
tau_study = prior(normal(0, 0.5), class = "sd") # tau_study
priors = list(
c(intercept100, SOO, NIV, tau_study),
c(intercept250, SOO, NIV, tau_study),
c(intercept500, SOO, NIV, tau_study),
c(intercept750, SOO, NIV, tau_study),
c(intercept1000, SOO, NIV, tau_study),
c(intercept2000, SOO, NIV, tau_study)
)
# Create list to store fits
# new_models = list()
#
# # Store the main model
#
# new_models[[1]] = m1
#
## Supressed because we saved the models already
#
## Run the loop, to fit other models with different
#
# for (i in 1:nrow(df)) {
#
# # Show what model is running
# print(i)
#
# # Skip first index because main model (m1) is already stored there
# k = i + 1
#
# # Re-fit the main model, but now changing the prior
# new_fit = update(m1,
# prior = priors[[i]],
#
# cores = parallel::detectCores(),
# sample_prior = T,
# chains = 4,
# control = list(adapt_delta = .99),
# backend = "cmdstanr",
# seed = 123)
#
# # Save models
# new_models[[k]] = new_fit
# }
# #
# saveRDS(new_models,
# here("03_updated_analyses_who", "output", "fits", "simulated_priors_models.rds"))
new_models = readRDS(here("03_updated_analyses_who",
"output",
"fits",
"simulated_priors_models.rds"))Diagnostic plots
# Run loop
for (i in 1:(nrow(df))) {
# Skip first index (main model)
k = i + 1
# Run custom function diag_plot, 1 per model
p = diag_plot(model = new_models[[k]],
pars_list = c("b_Intercept", "b_oxygenlow", "b_oxygenNIV", "sd_study__Intercept"),
ncol_trace = 4)
# Display plot
print(p)
# Add legend to show respective priors
lab = paste0(df %>% slice(i) %>% pull(SD) %>% round(2))
grid::grid.text(paste0("alpha_SD = ", lab), 0.25, .03, gp=grid::gpar(cex=1.3))
}
Posterior predictive check
# Run lopp
for (i in 1:(nrow(df))) {
# Skip first index (main model)
k = i + 1
p = pp_check(new_models[[k]], ndraws = 20)
print(p)
lab = paste0(df %>% slice(i) %>% pull(SD) %>% round(2))
grid::grid.text(paste0("alpha_SD = ", lab), .7, .9, gp=grid::gpar(cex=1.2))
}
Sensitivity analysis model
We will fit models assuming a priors for \(\alpha\) with the mean equal to \(log(1)\).
# Priors
intercept100 = prior_string(
paste0(
"normal(",
log(1), # Mean, 0
",",
df[[1,4]], # SD, total_arm = 100
")"),
class = "b", coef = "Intercept"
)
# b_Intercept ~ normal(0,0.28303369902181)
intercept250 = prior_string(
paste0(
"normal(",
log(1), # Mean, 0
",",
df[[2,4]], # SD, total_arm = 250
")"),
class = "b", coef = "Intercept")
# b_Intercept ~ normal(-0.261364764134408,0.179552672191073)
intercept500 = prior_string(paste0(
"normal(",
log(1), # Mean, 0
",",
df[[3,4]], # SD, total_arm = 500
")"),
class = "b", coef = "Intercept")
# b_Intercept ~ normal(-0.261364764134408,0.127092510741972)
intercept750 = prior_string(paste0(
"normal(",
log(1), # Mean, 0
",",
df[[4,4]], # SD, total_arm = 750
")"),
class = "b", coef = "Intercept")
# b_Intercept ~ normal(-0.261364764134408,0.103805945961058)
intercept1000 = prior_string(paste0(
"normal(",
log(1), # Mean, 0
",",
df[[5,4]], # SD, total_arm = 1000
")"),
class = "b", coef = "Intercept")
#b_Intercept ~ normal(-0.261364764134408,0.0899139032236548)
intercept2000 = prior_string(paste0(
"normal(",
log(1), # Mean, 0
",",
df[[6,4]], # SD, total_arm = 2000
")"),
class = "b", coef = "Intercept")
# b_Intercept ~ normal(-0.261364764134408,0.0635949873256457)
SOO = prior(normal(0, 1.5), class = "b", coef = "oxygenlow") # beta_Simple oxygen only
NIV = prior(normal(0, 1.5), class = "b", coef = "oxygenNIV") # beta_NIV
tau_study = prior(normal(0, 0.5), class = "sd") # tau_study
priors = list(
c(intercept100, SOO, NIV, tau_study),
c(intercept250, SOO, NIV, tau_study),
c(intercept500, SOO, NIV, tau_study),
c(intercept750, SOO, NIV, tau_study),
c(intercept1000, SOO, NIV, tau_study),
c(intercept2000, SOO, NIV, tau_study)
)
# Create list to store fits
# new_models = list()
# #
# # # Store the main model
# #
# new_models[[1]] = m1
# #
# # ## Supressed because we saved the models already
# #
# # # Run the loop, to fit other models with different
# #
# for (i in 1:nrow(df)) {
#
# # Show what model is running
# print(i)
#
# # Skip first index because main model (m1) is already stored there
# k = i + 1
#
# # Re-fit the main model, but now changing the prior
# new_fit = update(m1,
# prior = priors[[i]],
#
# cores = parallel::detectCores(),
# sample_prior = T,
# chains = 4,
# control = list(adapt_delta = .99),
# backend = "cmdstanr",
# seed = 123)
#
# # Save models
# new_models[[k]] = new_fit
# }
# #
# saveRDS(new_models,
# here("03_updated_analyses_who", "output", "fits",
# "sensitivity_simulated_priors_models.rds"))
new_models = readRDS(here("03_updated_analyses_who",
"output",
"fits",
"sensitivity_simulated_priors_models.rds"))Diagnostic plots
# Run loop
for (i in 1:(nrow(df))) {
# Skip first index (main model)
k = i + 1
# Run custom function diag_plot, 1 per model
p = diag_plot(model = new_models[[k]],
pars_list = c("b_Intercept", "b_oxygenlow", "b_oxygenNIV", "sd_study__Intercept"),
ncol_trace = 4)
# Display plot
print(p)
# Add legend to show respective priors
lab = paste0(df %>% slice(i) %>% pull(SD) %>% round(2))
grid::grid.text(paste0("alpha_SD = ", lab), 0.25, .03, gp=grid::gpar(cex=1.3))
}
Posterior predictive check
# Run lopp
for (i in 1:(nrow(df))) {
# Skip first index (main model)
k = i + 1
p = pp_check(new_models[[k]], ndraws = 20)
print(p)
lab = paste0(df %>% slice(i) %>% pull(SD) %>% round(2))
grid::grid.text(paste0("alpha_SD = ", lab), .7, .9, gp=grid::gpar(cex=1.2))
}
sessionInfo()## R version 4.0.4 (2021-02-15)
## Platform: x86_64-apple-darwin17.0 (64-bit)
## Running under: macOS Catalina 10.15.7
##
## Matrix products: default
## BLAS: /Library/Frameworks/R.framework/Versions/4.0/Resources/lib/libRblas.dylib
## LAPACK: /Library/Frameworks/R.framework/Versions/4.0/Resources/lib/libRlapack.dylib
##
## locale:
## [1] en_US.UTF-8/en_US.UTF-8/en_US.UTF-8/C/en_US.UTF-8/en_US.UTF-8
##
## attached base packages:
## [1] stats graphics grDevices datasets utils methods base
##
## other attached packages:
## [1] bayesplot_1.8.1 cowplot_1.1.1 patchwork_1.1.1 flextable_0.6.8
## [5] PNWColors_0.1.0 tidybayes_3.0.1 ggdist_3.0.0 rio_0.5.27
## [9] here_1.0.1 brms_2.16.1 Rcpp_1.0.7 gt_0.3.1
## [13] forcats_0.5.1 stringr_1.4.0 dplyr_1.0.7 purrr_0.3.4
## [17] readr_2.0.1 tidyr_1.1.3 tibble_3.1.4 ggplot2_3.3.5
## [21] tidyverse_1.3.1 pacman_0.5.1
##
## loaded via a namespace (and not attached):
## [1] uuid_0.1-4 readxl_1.3.1 backports_1.2.1
## [4] systemfonts_1.0.2 plyr_1.8.6 igraph_1.2.6
## [7] svUnit_1.0.6 splines_4.0.4 crosstalk_1.1.1
## [10] rstantools_2.1.1 inline_0.3.19 digest_0.6.27
## [13] htmltools_0.5.2 rsconnect_0.8.24 fansi_0.5.0
## [16] magrittr_2.0.1 checkmate_2.0.0 openxlsx_4.2.4
## [19] tzdb_0.1.2 modelr_0.1.8 RcppParallel_5.1.4
## [22] matrixStats_0.60.1 officer_0.4.0 xts_0.12.1
## [25] prettyunits_1.1.1 colorspace_2.0-2 rvest_1.0.1
## [28] haven_2.4.3 xfun_0.25 callr_3.7.0
## [31] crayon_1.4.1 jsonlite_1.7.2 lme4_1.1-27.1
## [34] zoo_1.8-9 glue_1.4.2 gtable_0.3.0
## [37] V8_3.4.2 distributional_0.2.2 pkgbuild_1.2.0
## [40] rstan_2.21.2 abind_1.4-5 scales_1.1.1
## [43] mvtnorm_1.1-2 DBI_1.1.1 miniUI_0.1.1.1
## [46] xtable_1.8-4 foreign_0.8-81 stats4_4.0.4
## [49] StanHeaders_2.21.0-7 DT_0.19 htmlwidgets_1.5.4
## [52] httr_1.4.2 threejs_0.3.3 arrayhelpers_1.1-0
## [55] posterior_1.1.0 ellipsis_0.3.2 pkgconfig_2.0.3
## [58] loo_2.4.1 farver_2.1.0 sass_0.4.0
## [61] dbplyr_2.1.1 utf8_1.2.2 labeling_0.4.2
## [64] tidyselect_1.1.1 rlang_0.4.11 reshape2_1.4.4
## [67] later_1.3.0 munsell_0.5.0 cellranger_1.1.0
## [70] tools_4.0.4 cli_3.0.1 generics_0.1.0
## [73] broom_0.7.9 ggridges_0.5.3 evaluate_0.14
## [76] fastmap_1.1.0 yaml_2.2.1 processx_3.5.2
## [79] knitr_1.34 fs_1.5.0 zip_2.2.0
## [82] nlme_3.1-152 mime_0.11 projpred_2.0.2
## [85] xml2_1.3.2 compiler_4.0.4 shinythemes_1.2.0
## [88] rstudioapi_0.13 gamm4_0.2-6 curl_4.3.2
## [91] reprex_2.0.1 bslib_0.3.0 stringi_1.7.4
## [94] highr_0.9 ps_1.6.0 Brobdingnag_1.2-6
## [97] gdtools_0.2.3 lattice_0.20-44 Matrix_1.3-4
## [100] nloptr_1.2.2.2 markdown_1.1 shinyjs_2.0.0
## [103] tensorA_0.36.2 vctrs_0.3.8 pillar_1.6.2
## [106] lifecycle_1.0.0 jquerylib_0.1.4 bridgesampling_1.1-2
## [109] data.table_1.14.0 httpuv_1.6.3 R6_2.5.1
## [112] promises_1.2.0.1 renv_0.14.0 gridExtra_2.3
## [115] codetools_0.2-18 boot_1.3-28 colourpicker_1.1.0
## [118] MASS_7.3-54 gtools_3.9.2 assertthat_0.2.1
## [121] rprojroot_2.0.2 withr_2.4.2 shinystan_2.5.0
## [124] mgcv_1.8-36 parallel_4.0.4 hms_1.1.0
## [127] grid_4.0.4 coda_0.19-4 minqa_1.2.4
## [130] cmdstanr_0.4.0 rmarkdown_2.10 shiny_1.6.0
## [133] lubridate_1.7.10 base64enc_0.1-3 dygraphs_1.1.1.6